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In the harmonic description of general relativity, the principle part of Einstein's equations reduces 
to 10 curved space wave equations for the components of the space-time metric. We present theorems 
regarding the stability of several evolution-boundary algorithms for such equations when treated in 
second order differential form. The theorems apply to a model black hole space-time consisting 
of a spacelike inner boundary excising the singularity, a timelike outer boundary and a horizon 
in between. These algorithms are implemented as stable, convergent numerical codes and their 
performance is compared in a 2-dimensional excision problem. 
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INTRODUCTION 



A primary goal of numerical relativity is the computation of gravitational radiation waveforms from binary black 
holes. Radiation produced in the inspiral and merger of binary black holes is expected to provide a strong signal for 
gravitational wave observatories. However, the simulation of black holes has proved to be a difficult computational 
problem. The importance of this challenging problem has recently spurred a fertile interaction between numerical 
relativity and computational mathematics. The classic computational treatment of hyperbolic systems has been 
directed at fluid dynamics and has been based upon first differential order systems. Certain formulations of Einstein's 
equations take a more natural second order form, notably the harmonic formulation |l],l2] for which well-posedness of 
the Cauchy problem was first established |j|. Here we present theorems regarding the stability of several evolution- 
boundary algorithms for such second order systems which have direct application to the black hole problem. 

Harmonic coordinates x a — (t,x l ) — (t,x,y,z) have only recently been used in designing numerical codes |j,L|ui 
^S. ' Q) l3> t» llfl fill Il2 | . They satisfy the curved space wave equation 

U . ng3 M := ' da[ ^- gg ^ dpx ^ = o. (i.i) 
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(3JT), In harmonic coordinates, Einstein's equations reduce to 10 quasilinear wave equations for the components of the 

metric, 



^ : u gg r = s^, (i.2) 

H 

where S^ u are nonlinear terms which do not enter the principle part. Thus the scalar wave equation 

g ap d a d l3 u = 0, (1.3) 

which has the same principle part, provides a fundamental testing ground for designing algorithms to treat the 
nonlinear gravitational problem IJ1.2J1 . In a previous study J.13J, we used this scalar equation to develop evolution and 
boundary algorithms for a model one dimensional black hole excision problem. Here we extend these results to two 
dimensions. While the extension to 2D involves substantial new features, the generalization from 2D to 3D is quite 
straightforward. Thus our results are immediately applicable to algorithms for the harmonic gravitational equations 
(ll.2|l . as well as their generalization to include harmonic gauge forcing terms [Tj| and other related generalizations 
such as the ZA formulation pj|. 

We treat 111. .'it in the second order differential form, which has advantages for both computational efficiency and 
accuracy over first order formulations PJJUJl- Although the system can be reduced to first order symmetric hyperbolic 
form Jig, this has the disadvantage of introducing auxiliary variables with their associated constraints and boundary 
conditions. The second order form is also best suited to the analogous wave equations governing elasticity and 
acoustics. Elasticity theory is governed by a coupled system of wave equations which for simple cases is similar to 
111.3(1 . in which the spatial components g %3 are determined by the elastic moduli. In fact, some of the techniques utilized 
here have been developed in a recent computational study of the wave equations governing an elastic body []jj]. The 
new ingredient introduced in the wave equation 1I1.3JI arises from the non-vanishing mixed space-time derivatives 
arising from the components g lt . Such terms do not ordinarily appear in the wave equations governing elasticity 



theory because they are treated in the rest frame of the body but they would necessarily arise in treating acoustic 
waves propagating in a medium with nonuniform macroscopic motion. In general relativity, these mixed space-time 
components of the metric correspond to a non- vanishing "shift", which is an essential feature of the black hole problem. 
In the standard 3 + 1 description of space-time J2CJ], the Cauchy hypersurfaces t = constant are required to be spacelike 
so that they have a length element with Euclidean signature 

d£ 2 =h ij dx i dx j . (1.4) 

The inverse spatial metric h % i , satisfying h % ^h^ = 5 l k , is related to the spatial components of the 4- metric determining 
the wave operator by 

h ij =g ij --^p i j3 j , (1.5) 

g 

where f3 l are the components of the shift. Here, the spacelike character of the Cauchy hypersurfaces requires that 
g tt < 0. 

The wave equation with shift has not received a great deal of attention outside of recent work in general relativity 0, 
Iirlll^.l2lll22.l23|. Even before the computational problem is attempted, new mathematical features introduced by 
the shift must be incorporated into the formulation of a well-posed initial-boundary value problem. The operator 
h^didj is by construction an elliptic operator defined by the spatial metric of the Cauchy hypersurfaces. However, 
the operator g^didj is elliptic only when the shift is sufficiently small. The elliptic case arises when the operator dt 
is timelike, i.e. when the evolution proceeds in a timelike (subluminal) spacetime direction. 

Without loss of generality, we set g tt = — 1 and write the 2D version of (II, 3|) as 

d 2 t - 2{p*d a + f3 y d y )d t - K - (3 x p x )dl - (ci - p»p y )d 2 y - 2(h - f3 x py)d x d y \ u = o, (1.6) 

where h xx = oi, h xy = b\ and h vv = c\ . The Euclidean property of h^ requires 

ai>0, ci>0, aici-bf>0. (1.7) 

The components of g % i are g xx — a = a\ — /3 x f3 x , g vv = c = c\ — (3 y (i v and g xy = b = b\ — [3 X fi v . In the subluminal 
case when g^didj is an elliptic operator, the simplest second order accurate difference approximation to Hl.fil) is 



W := I df - 2(P X D QX + P»D 0y )d t - aD +x D_ x - cD +v D_ v - 2bD 0x D 0y \u = 0. (1.8) 

(Here Dqi, D + i and Z?_i are, respectively, the centered, forward and backward difference operators in the ^-direction 
defined in Sec. lIIIJl . This leads to stable evolution-boundary algorithms for Dirichlet, Neumann, Sommerfeld or other 
dissipative boundary conditions. Stability was established for the ID case using a semi-discrete energy norm in [13j . 
and this was generalized using the discrete energy method to the full 3D case in (lfl Il2l | . 

The W-algorithm i|1.8|l is unstable when the shift is sufficiently large so that g l% < (for any diagonal component). 
This occurs in one of the strategies for avoiding problems with the singularity which ultimately forms inside a black 
hole. In this strategy, the singularity is "excised" by surrounding it with a spacelike inner boundary. The evolution 
direction which is adopted to this spacelike boundary is superluminal, so that g v no longer has Euclidean signature. In 
Sec. IIIII we establish the stability of several different second order, evolution-boundary algorithms for this superluminal 
case. For the Cauchy problem, we establish stability for a general system of wave equations in s spatial dimensions 
so that the results may be immediately applied to other second order systems such as elasticity theory or acoustics. 
For the boundary, we specialize our treatment to scalar equations in 2D in order to simplify the notation, but the 
extension to general systems in sD is straightforward. In particular, our results apply directly to the 3D harmonic 
evolution of black holes. 

The analysis of the initial-boundary problem for ijl,6H in Sec. IIIII makes evident that the above geometric properties 
of the wave equation have a mathematical analogue which results independently from a consideration of the well- 
posedness of the problem. The geometrical and analytical approaches are complementary and provide a good meeting 
ground for the ideas of numerical relativity and computational mathematics. While the main concern of numerical 
relativity is the black hole problem, the stability theorems for the finite difference algorithms developed for the model 
problems considered here provide a firm basis for attacking this problem with the harmonic Einstein system 111, 211 . 

In Sec. lIVI we compare the performance of the algorithms for the superluminal case in a problem without boundaries. 
In Sec.0 we simulate a simple 2D model of the excision problem in which the inner boundary S is spacelike and the 
outer boundary T is timelike. Between the boundaries the operator g^didj goes from non-elliptic to elliptic along a 



curve TC where det(g 1 ^) = 0. The metric is chosen so that no characteristics can leave the inner region between S and 
H, so that H mimics the role of a horizon. The global simulation of l|1.6|l in the region bounded by T and S is achieved 
with a blended evolution algorithm. A stable superluminal algorithm is used in an inner region between S and H. 
In the exterior region, this superluminal algorithm is blended to the W^-algorithm IJ1.8J1 . so that the W^-algorithm is 
used to treat the outer boundary T. This model excision problem involves many of the mathematical difficulties in 
the full gravitational case. We begin in Sec. [H] with some simple examples which illustrate the problem, its potential 
pitfalls and how to avoid them. 

II. SOME SUBTLETIES ASSOCIATED WITH THE WAVE EQUATION WITH SHIFT 

In an inertial coordinate system x a = (£, x 1 ) (in units where the velocity of light c = 1), the wave equation which 
governs special relativistic physics, 

fl£-<SfcW)u = 0, (2.1) 



does not contain a shift. The invariance of the velocity of light results from the property that this wave equation 
retains the same form 12.1JI under a Lorentz transformation, 

t' = , 1 (i-5 iJ i &) 

-AS* -pi) 



fi 2 = 8 kl (3 k f3 l , (2.2) 

to another inertial coordinate system with relative motion. In this way special relativity resolves the dilemma with 
experiment that under a Galilean transformation, 

t = i 
x l = (.t'-/3 4 £), (2.3) 

(12, 1|) gives rise to the shifted wave equation 

d 2 t - 2(3% - (6 i: > - I3 l (3 1 )d l d^\u = (2.4) 

whose solutions propagate with coordinate speeds in the range |1 ± (3\ (where (3 2 = S^^ft). This raises the question: 
why does the wave equation with shift arise in general relativity? 

In fact, although there are no preferred inertial coordinates in general relativity, in any sufficiently small spacetime 
region it is always possible to introduce Gaussian coordinates in which the wave equation 111. -Ill reduces to the shift-free 
form 

(<9 t 2 - h ij didj)u = 0. (2.5) 

The problem here is that in Gaussian coordinates the worldlines x r = constant are geodesies, i.e. the worldlines of 
freely falling observers, which can be focused by the attractive nature of gravity to produce coordinate singularities. 
This can occur on a short time scale in a strong gravitational field. 

Another reason for introducing a shift is the simplicity of harmonic coordinates in reducing Einstein's equations 
into the hyperbolic form ljl.211 . Since the shift components g lt satisfy a coupled system of nonlinear wave equations, 
even if they were initialized with vanishing Cauchy data they would in general evolve to be non-zero. This cannot 
be avoided by introducing a harmonic gauge forcing term, of the form Ux a — F a , without choosing the forcing term 
F a to depend upon the derivatives of the metric d^g^ . This in turn jeopardizes the hyperbolic form of the reduced 
Einstein equations and the well-posedness of the Cauchy problem fl4| . 

Another reason for introducing a shift arises in the simulation of black holes. Once a black hole of mass M has 
formed there is at most a proper time of order M (in gravitational units) until a physical singularity is encountered. 
On the other hand, a simulation which provides gravitational waveforms of physical interest typically requires an 
evolution for a proper time of more than 100A/ in the exterior region. One strategy for accomplishing this is to excise 



the singularity by surrounding it with a spacelike inner boundary for the simulation domain, i.e. an inner boundary 
which moves at superluminal speed. If the evolution tracks the inner boundary then a superluminal shift must be 
used. 

This can be illustrated by a spherically symmetric Schwarzschild black hole for which the wave equation ll.fljl 
becomes 

(1 + 2 J£ )d ? _ ^L dtdr _ (1 _ 2 J± )d * r 1 {d l + -±-dl))u = 0, (2.6) 

r r r r* sin 6 v J 

in ingoing Eddington-Finkelstein coordinates. Here the evolution takes place on the spacelike Cauchy hypersurfaces 
t = const which are non-singular for r > 0. The black hole is located at r = 2M, which is a characteristic hypersurface 
with the horizon property that no characteristics leave the region r < 2M. The singularity is excised by evolving in 
a domain R\ < r < R 2 , where < R\ < 2M and Ri » 2M. The shift has the radial component 

(3 r = — i^- > 0. (2.7) 

1 + 1M 

The change in sign of the coefficient of d^ in passing inside the horizon does not change the hyperbolicity of the wave 
equation but it changes its mathematical properties. Outside the horizon, the curves of constant (r, 9, <j>) are timelike, 
as well as the outer boundary r = R 2 . In the outer region 2M < r < fc , the VF-algorithm 1|1.8I) provides a stable 
second order evolution-boundary algorithm for the wave equation [jjj, HH llil • 

Inside the horizon, the t-direction, as well as the inner boundary r = Ri is spacelike, i.e. evolution on a grid with 
constant (r, 6, <f>) proceeds outside the light cone. This effects the mathematical properties of the wave equation. As 
a result, in this domain, the VF-algorithm is unstable. The alternative algorithms presented in Sec. IIII J ar e stable 
inside the horizon. But the ^-algorithm has better accuracy than these algorithms in the exterior region [lj|. In the 
simulation of the model excision problem in Sec. a stable algorithm for the superluminal regime is blended to the 
Vy-algorithm in the exterior. 

The Schwarzschild horizon has the property that characteristics can not exit from inside, but can enter from the 
outside. Near the horizon, the radial part of Schwarzschild wave equation Ij2.fi 11 has the same qualitative features as 
the wave equation 

(d t -d x )(d t + xd x )u = 0, (2.8) 



which has a horizon x = 0. One set of characteristics of 12.8JI cross the horizon at x = in the negative x-direction. 
The other set of characteristics are tangent to the horizon and diverge away on either side. An observer at x > 
cannot see beyond the horizon at x = 0. This is the situation which we dealt with in a model ID excision problem (l3j 
whose treatment we generalize to 2D in Sec. However, it should be emphasized that the related equation 

(d t : - d x ){dt - xd x )u = (2.9) 

has a different mathematical character. Although \2M\ is also hyperbolic and has a well-posed Cauchy problem, one 
set of characteristics converge toward the horizon at x — 0. These characteristics approach each other exponentially 
fast and, in general, the gradients become exponentially large near x = 0. This would lead to the focusing of a wave 
into formation of a shock. Although we do not treat this case in this paper, it is important to bear in mind that it 
would require different methods. 

Boundaries introduce additional subtleties. First consider a timelike boundary, similar to the outer boundary 
r = R-2 > 2M for the Schwarzschild wave equation H2.fijl . Since the evolution is timelike in the neighborhood of the 
boundary, the M^-algorithm can be used. The stability of dissipative boundary conditions for the VK-algorithm was 
established for ID in [JJ| and extended to 3D in [12[ by means of a semi-discrete energy method. However, such 
an energy estimate does not preclude exponential growth of a wave traveling between two boundaries. A simple 
example Q arises from the repetitive blue shifting of a wave packet in special relativity reflecting back and forth 
between two plane boundaries, whose velocities ±i> are controlled to be always toward the packet during reflection. 
After many reflections the wave packet shrinks in size and its energy grows by a factor e 4aT , where T is measured 
in units of the crossing-time between reflections and v = tanh a. Dissipation must be used to control such growth of 
short wavelength error. 

It is instructive to interpret the boundary conditions on a wave in special relativity in the shifted coordinate system 
(12. 3|) where the boundary has fixed location but moves relative to the t = const Cauchy hypersurfaces. In the ID 
case, this gives rise to the half-plane problem 



dt-2pd x dt-(l-F)%)u = 0, (2.10) 



in the region x < (where we now write (3 X = (3) . There are two different frames in which the energy of the wave can 
be considered - the rest frame of the boundary and the rest frame intrinsic to the Cauchy hypersurfaces. In the rest 
frame of the boundary, the energy is 

E =\J dx((d t uf + (I - (3 2 )(d x u) 2 \ (2.11) 

and satisfies 

d t E = d t u ({I - I3 2 )d x + 0d?\ u|x=o. (2.12) 

In the case 1 < 1 , this energy provides a norm and the semi-discrete version of the flux-conservation law ( 12.1211 
provides the basis for establishing stable evolution-boundary algorithms for the VF-algorithm ijl,8Jl , Note the sign 
of P is important here in formulating a stable Neumann boundary condition. A homogeneous Neumann boundary 
condition takes the dissipative form 



(l-/3 A )d x +pd t \u = 0. (2.13) 

The familiar form d x u = implies dtE < and thus guarantees a well-posed problem only when /? > 0, i.e. only 
when the motion of the boundary is outward relative to the Cauchy hypersurfaces. 
The energy intrinsic to the Cauchy hypersurfaces, 



p 

2J-oo~~V 

provides a norm even in the superluminal case when j3 2 > 1. It satisfies 

P, a .. aa .,2 , P fa ^2 



Eo = l [ dx ((d t u - I3d x uf + (d x u) 2 ) , (2.14) 



d t E = {^(d t u-/3d x uy + ^(d xU y + (dtU-pd x u)d x u)\ x =o. (2.15) 

Thus, in the absence of a boundary, (I2.15|l would reduce to d t E = so that the Cauchy problem is well-posed for 
any j3. The energy analogous to -Bo is used in Sec. IIIII to establish well-posedness of the Cauchy problem and the 
stability of superluminal algorithms in the general multi-dimensional case. 

When < — 1, i.e. when the motion of the boundary is superluminal and directed toward the Cauchy hypersurfaces, 
it is easy to verify that J2.150 implies dtEo < so that there is always a loss of energy through the boundary. This 
is the case of a spacelike boundary through which all the characteristics leave, i.e. a pure "outflow" boundary. 
Stable algorithms for such a boundary are also given in Sec. IIIII for the higher dimensional case. Note that for 
P > 1 the boundary is also spacelike but now IJ2.15J) implies dtEo > 0. This is the pure "inflow" case, in which all the 
characteristics enter the boundary. This should not be considered in the context of an initial-boundary value problem, 
but as a pure Cauchy problem where the boundary represents a non-smooth extension of the Cauchy hypersurface. 

Further subtleties arise in treating co-orbiting, binary black holes. One strategy for the binary problem is to use a 
rotating coordinate system which co-orbits with the black holes. In the Schwarzschild case, the use of a coordinate 
<p = cf> — u>t rotating with angular velocity uu transforms the wave equation H2.6J1 into 

(1 + ™L m + ud f _ ^£ dtdr _ (i _ 2M )Q a _ 2_ { g 2 + 1&)) U = o. (2.16) 

r r r H sin 6 J 

Now the t-direction becomes spacelike in the region 

(l + ^.) r Vsin 2 0>l, (2.17) 

r 

which intersects the outer boundary r — R2 if R2 is sufficiently large. In that case, although the boundary remains 
timelike the evolution is superluminal so that the ^-algorithm is no longer stable. A stable algorithm for such a 
boundary problem has been established in the ID case [2j|. We will not consider the 2D version of this problem here. 
A common strategy for treating the binary black hole problem is to use a grid based upon Cartesian coordinates. 
This poses a problem in dealing with inner and outer boundaries with the spherical shapes natural to the problem. In 
other second order wave problems, such curved boundaries have been successfully treated by the embedded boun dary 
method [l3,l2J|- Another approach being explored in general relativity is to use multi-block grids |25l I26L I27L l2fj| . 
This is another problem which we defer to future work and do not consider here. 



III. ALGORITHMS FOR THE 2D SUPERLUMINAL PROBLEM 

In this section, we study a class of second order hyperbolic systems with shift which we will use in Sec. ^ to 
construct stable algorithms for a model 2D black hole excision problem. The excision problem is a strip problem 
with spacelike and timelike boundaries and a horizon in between. In the region where the shift is superluminal, the 
boundary is spacelike and where the shift is subluminal, the boundary is timelike. We replace this problem by Cauchy 
and half-space problems. The strip problem is well-posed if the corresponding Cauchy and half-space problems are 
well-posed |3(l |. 

For the Cauchy problem, we consider general systems of equations in s space dimensions to demonstrate that the 
results have applicability beyond numerical relativity. For the half-space problems, we only consider scalar equations in 
2D to simplify the notation. The generalization from scalar equations in 2D to systems in sD is quite straightforward. 

Here, we consider systems with constant coefficients. Systems with variable coefficients can be reduced to systems 
with constant coefficients by freezing the coefficients at all points. The problem with variable coefficients is strongly 
well-posed if the Kreiss condition holds uniformly for all problems with constant coefficients j^jj . 

In order to analyze and establish stable approximations we use the method of lines and reduce the system of 
partial differential equations to a system of ordinary differential equations in time on a spatial grid. We then apply 
two standard techniques: the energy method and mode analysis. The stability of the semi-discrete approximation 
implies the stability of the totally discretized method for most standard methods of lines J22|, e.g. with the use of a 
Runge-Kutta time integrator. 

A. The Cauchy Problem 

We consider the Cauchy problem for a second order system with constant (possibly complex) coefficients in s space 
dimensions, 



u tt 



with the initial conditions 



^A jk ——n:=P {d/dx)vL, x = {x x , ...,x 8 ) e R s , t > 0, (3.1) 

j,k=l 3 



u(x,t = 0) = f(x), ut(x,t = 0)=g(x), u,f,gG<C". (3.2) 



(We abbreviate d a u — u a where confusion does not arise.) Here, for each (J, k), A^ are constant Hermitian matrices 
G C n '™, and the data f = f(x) and g = g(x) are smooth and 1-periodic in each Xj, j = 1, ...,s. The solution 
u = u(x, t) is then smooth and 1-periodic in each Xj. Moreover, we consider solutions with J Rs udx = 0. 

We assume that the Hermitian operator Pq in (13. 1|) is elliptic, i. e. there exists a positive constant S such that 

s 

J2 A ^£u > m\ 2 i (3.3) 

j,k=i 

for all vectors (eR s . Here / is the n x n identity matrix. 
We introduce a shift by 

x = x + /3t, x=(x 1 ,...,x s )GM s , /3=(/3\...,/? s )gR s , /3 j >0, 

and obtain the shifted system 

u tt = 2Pi (d/dx)u t - Pi (d/dx)u + P (d/dx)u. (3.4) 

Here Pi is a scalar operator, 

Pi(d/dx) = J2/3^. 

3=1 ^ 

Theorem 1. The Cauchy problem for i'i.H is well-posed. 



Proof. If we set v = u t — P\(d/dx)u, we get the first order system 



We Fourier transform 1'i.ot and get 

J),- am (J) + (-a mO(J)- """'■ (M) 

where u = (u>i, . . . , lo s ) £ Hl s and P\(iuj) = i Yfj=i J Uj and Pq(u) = 5ZJ &=i AjkWjUlk- Since Po is elliptic, we have 
Pa = Pq > 5o\u>\ 2 I, for some 5q > 0. We can then introduce new variables 





/ 


'u\ 




I 


w 


= T( 


,♦) 


, r = 


5 -1/2 

I Po 



(3.7) 
and, since T and Pi commute, we obtain from (13.611 

w t = A(iw)w + ( A J/o P ° (W) | w := Sw. (3.8) 

V -ft » o ; 

Since the matrix S is skew Hermitian, S* = —S, we obtain 

4-||w|| 2 = (5w,w) + Mr, Sir) = (w, (S* + 5)w) = 0. 
at 

Therefore, by Parseval's relation, there is an energy estimate and the Cauchy problem is well-posed J21|. □ 

Now we show how to construct stable finite difference approximations to (13, 4|) , We leave time continuous and use 
the method of lines. For brevity, we treat the case /J- 7 > 0. 

Let hj = 1/Nj, j = l,...,s, denote spatial gridlengths, where Nj are natural numbers. For any multi-index 
v = (vi,...,v B ) £ Z s , let Xy = (h\v\, . . . ,h 8 v 8 ) denote the corresponding gridpoint. We consider gridfunctions 
u„ := u,/(Xy, t) approximating u(x„, t) and introduce a translation operator Ej in the j-th coordinate by 

Pju„ = u^x,, + phjej,t), peZ, 

where ej = (0, . . . , 0, 1, 0, . . . , 0) is the vector containing a 1 in the j-th position and zeros elsewhere. We then define 
the forward, backward, and the central difference operators in the j-th coordinate direction by 

hjD +j = E}-E°, h j D_ j = E G j -Ej\ 2D 0j =D +J +D^. 

We approximate (|3.4|l by 

Uutt = 2 Pl (D)u ut - pl(D)u u + Po (D)u v , (3.9) 

where po (D) is the centered approximation 

s s 

p {D) = ^A jj D +j D_ j + Y, A jk D 0j D k, (3.10) 

3 = 1 j^k=l 

and Pi(D) is any one of the following approximations: 
1) Centered approximation, 

s 

Pl (D)=J2^D 0j , (3.11) 



2) First order accurate one-sided approximation, 



Pl {D) = Y,P j D +h (3.12) 



3) Second order accurate one-sided approximation, 



Pl {D)=J20 J D pj , (3.13) 

3=1 



where 

n . — n . . 

2 +•* 



D pj = D +j - ^D* (3.14) 



Remark. It is not necessary to assume that p J > in l|3.11ll - l|3.13ll . In general, we can use - — ^— l D + j + - — ^—^-D-j in 

13.12JI . For the second order one-sided approximation (I3.13|l . we replace D + j and D_j by D p j and D m j = D-j + -j-D 2 _ :j , 
respectively. 



Theorem 2. The approximation IJ3.9H is stable. 

Proof. As in the continuum case, we write 113. 9JI as a first order system and Fourier transform to get 



where 



° 4 £ \ 

Po = Yl A n ; T2 sin2 ~2 + Yl Ajkj-j-smtjSmZk, ^=Ujhj, |&| < tt, (3.16) 



i=i i#fe=i 

and pi is one of the following 



p 1 =J2j3 j ±-ism£ j , (3.17) 



Pi 



j=i 3 



Pi = ^/3 J 'iusin^-2sin 2 |) 



^i. ^ sin ^(l - 2sin 2 |) - 4sin 4 |) , (3.19) 



corresponding to H3.llll - H3.13l) . 
Since 



we have 



£ £ £ £ 

sin 2 £ = 4 sin 2 - cos 2 - = 4 sin 2 - - 4 sin 4 - , (3.20) 

s 2 2 2 2 y ' 



Po = J2 A * hk~ Shl ^ Sin & + £ AjJ /? Shl4 2' (3 - 21) 

From the ellipticity condition l|3.3|l it follows that po is positive. As £j — ► we have sm^j/hj — > Wj. Therefore, the 
first sum in IJ3.21J) is strictly positive. When |£j| = 7r, the first sum in IJ3.21I) is zero but the second sum is not because 
sin t=J- 7^ 0. Therefore po is positive definite, and we can use the same transformation as in IJ3.7JI and write l|3.15Jl as 



o f>y\uj) 

-po /2 M o 



w t =p!W+( „ i; u 2 N ^ V ^ jw. (3.22) 



The second term on the right hand side of (I3,22|l is again skew Hermitian and has no influence on the stability. 
Thus, we need only consider 

which consists of difference approximations of scalar equations of the above type. To show that the approximations 
(I3,llj) - (l3,13j) are stable, we set u = e xt u and get A = pi. By Il3.17jl - j3.19jl . we have 5RA < and there are no 
exponentially growing modes. □ 

The approximation IJ3.9H involves a wide stencil. Therefore extra boundary conditions (ghost points) are required 
and the resulting accuracy is less than with a more compact stencil. In order to investigate other approximations 
with a more compact stencil, we write (13. 411 as 

u tt = 2P 1 (d/dx)ut + P(d/dx)u, P(d/dx) = P (d/dx) - Pl{d/dx) (3.23) 

and approximate it by 

u vtt = 2p 1 (D)u vt +p(D)u v , (3.24) 

where p\ (D) is given by (|3.11|l and p(D) is the centered approximation 

s s 

p(D) = Y,{ A u-P j2 ) D +i D -J+ E (A jk -^0 k )D Oj D ok . (3.25) 

Theorem 3. The approximation 13.2411 is stable if Ajj — ft > 0. 
Proof. We write Q3.24H as 

u 1/tt = 2pi{D)u 1/t -p 2 1 (D)u 1/ + q{D)u v , q{D) = p{D) + p\{D) . (3.26) 

We use the relation D + jD-j = D^j — ^-D+jD^_j and write 



1(D) = ]T A jk D Qj D ok - -Y.{A n -ft 2 )h]D\ 3 Dl r 



j,k=i j=i 

In the same way as in the continuum case, we write (I3.2fijl as a first order system and Fourier transform to get 

where 



v J \v/ y — q J \ v 



E A ^-Th~ Sin ^' Sin?fe + E ( A m - ^ 2 ) ^ sin4 9 ' & = W A> &l ^ "■ ( 3 - 28 ) 



9 ^ " jK h,h k - "^ ' ^ \" jj ) m - 

3 ,k=l J i=l J 



2 



and pi is given by (I3.17J1 . By the ellipticity condition 13. .'it , it is clear that q is a positive definite matrix if Ajj — ft > 
Therefore, we can use the same transformation as in IJ3.7II and write IJ3.27J1 as 

o <f/ 2 H 

-q 1 / 2 ^) 



w t =prW+( _m /2 , , y n V ; )w. (3.29) 



The second term on the right hand side of (I3.29|i is again skew Hermitian and has no influence on the stability. 
Thus, we need only to consider 

w t =piw, 

and the stability follows in the same way as in Theorem 2. □ 
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Remark. If the operator P is elliptic, we have Ajj — ft > 0, and by Theorem 3 the approximation 1.1.2411 is stable. 

However, it is possible to have Ajj — ft > while P is non-elliptic. In this case, the approximation l.'i.24t remains 
stable when P is non-elliptic. In other words, the stability of i't.2411 does not depend upon the coefficients of mixed 
derivatives Aj k , j ^ k. 

Remark. In the scalar case, (|H.24|) reduces to the l^-algorithm (II, 8|) , 

In the excision problem, we use the subluminal algorithm IJ3.24J) in the subluminal region where Ajj — ft > 0. In 
the superluminal region where the shift ft is large so that Ajj —ft < 0, we use the superluminal algorithm 13.9JI 
instead. We need then a prescription for switching from one algorithm to the other. There are two distinct ways to 
do this. One is to make a sharp switch between the algorithms where the transition from superluminal to subluminal 
region takes place. The other, used in []j|], is to introduce a smooth, monotonic blending function and use a blended 
algorithm, which turns into the superluminal algorithm inside the superluminal region and reduces monotonically to 
the subluminal algorithm in the outside. For this purpose, note that the superluminal algorithm remains stable in 
the subluminal region. 

As a further alternative to the above approximations, we can approximate l|3.23|l by adding a fourth differential 
order term 

u vU = 2p 1 (D)u lJt + p(D)u, y - Q(£>H, (3.30) 

where p\ (D) is given by (|3.11|l and 

i s 
Q(D) = jYl a i h2 3 D +3 D2 -v a J ^ °- ( 3 - 31 ) 

The motivation for adding such a fourth order term is to modify the matrix q in Ij3.28|) so that it becomes positive 

2 1 

definite even if Ajj —ft < 0. When Ajj —ft" > 0, the matrix q is positive definite and this added term is 
unnecessary. We can take advantage of this by embedding the switch or blending function in the choice of ctj , with 
olj = in the outer region. 

Theorem 4. The approximation H3.30|) is stable if Ajj + ajl > ft I. 
Proof. We use the relation Dqj = D + jD-j + -j-D+jD^j and write 

S 1 s 

p(D) = £ (A jk ftft)D 0l D 0k - J^(Ajj - ft 2 )h 2 jD 2 +j D 2 _j 

j,k=i j=i 



= -p\{D) + J2 A 3k D 0] D 0k \ J2(Ajj - ft^hpl^. 

3 ,k=l 3=1 

We can then write 113. 3011 as 

u vtt = 2 Pl {D)u ut - pi(D)u u + q{D)u v , (3.32) 

where 

q(D) = J2 A jk D 0j D ok - - J2(A n - ft 2 + aj)h 2 DljD 2 _j. (3.33) 

j,k=l 3 = 1 

In the same way as before, we write H3.32J1 as a first order system and Fourier transform to get 

.:),-*C) + (4,i)(:)- 

where 

q = Y, A Jk J^~ sin ^ sin & + ^ ( A oi ~ F ^ + a 3 ^2 sin4 f • 
j,k=i ° " j=i J 

If Ajj + ajl > ft I then, because of ellipticity, q is positive definite and stability follows in the same way as before. 
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B. Half-plane Problems 

We consider the scalar wave equation with constant coefficients in two space dimensions, 

u~ t i = oi u xx + 2&i u xy + c\ Uyy •— P u. (3.35) 

In the moving coordinate system, t = t, x = x — (3 x t, y = y — (3 v t, with f} x ,f] v > 0, we get the shifted wave equation, 

u tt = ^{P x u x t + /3 v u yt ) + au xx + 2bu xy + cu yy := 2P 1 u t + Pu. (3.36) 

Here the coefficients a — a± — (3 x2 , b = b x — f3 x (3 y , and c = C\ — [3 x2 are assumed to be constant. Moreover, we assume 
that the space operator Pq in l|3.35|l is elliptic, namely a\ > and c\ > and b\ < a\C\. Therefore, by Theorem 1, 
the Cauchy problem for Il3.36p is well-posed. 
We consider Ij3.36ll in the half-space 

< x < oo, — oo < y < oo, t > 

and we assume that u is 1-periodic in y. The number of boundary conditions needed at x = is equal to the number 
of outgoing characteristics of the equation u u = 2f3 x u xt + au xx . We consider two distinct half-plane problems 
determined by the coefficients of the operator P. 

Half-plane problem I: If a > and b 2 < ac, then the operator P is elliptic and one boundary condition is needed 
at x = 0. In the excision problem, this is the case of subluminal shift with a timelike boundary. 

Half-plane problem II: If a < 0, then the operator P is non-elliptic. In the excision problem, this is the 
case of a superluminal shift with a spacelike boundary. 

1. Half-plane Problem I (Subluminal Case) 

This is the problem treated in [l2| by the energy method. In the present context of IJ3.36II . the energy is given by 

E = |K|| 2 + all^f + 2b(u x ,u v ) + c\\u y \\ 2 (3.37) 

in terms of the L2 scalar product and the corresponding norm 

(v,w) = / vwdxdy, \\v\\ 2 = (v, v). (3.38) 

JQ JO 

If u solves 13.36JI , then integration by parts gives 

d t E= -2ut{[3 x ut + au x + buy) . (3.39) 

x— 

Any boundary condition satisfying the dissipative condition dtE < gives an energy estimate sufficient to establish 
the well-posedness of the Cauchy problem, including the Dirichlet condition 

u t (0,y,t) = (3.40) 

and the Neumann condition 

[3 x u t (0, y, t) + au x (0, y, t) + bu y (0, y,t) = (3.41) 

for which energy is conserved. 

As difference approximation for the half-plane problem, we use l|3.24|i . which in the present case reduces to the 
ly-algorithm IJ1.8J1 . By introducing a discrete energy norm and using summation by parts, a discrete version of 113. 39J) 
has been used to establish stability of the finite difference problem. For details we refer to |12 |. 
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2. Half-plane Problem II (Superluminal Case) 

To investigate the well-posedness of the continuum problem, we use mode analysis. We apply a Laplace transfor- 
mation in t and Fourier transformation in y. 

Theorem 5. The half-plane problem l|3.36l) with a < is well-posed. 
Proof. By substituting u = u(x) e st+iujy , sGC^el, into (13.3611 we obtain 

au xx + (2iboj + 23 x s)u x + (2i/3 v ujs - s 2 - cuj 2 )u = 0. (3.42) 

The general solution to the ordinary differential equation <3.42ll is of the form u{x) = o\e KlX + a2e K2X , where K\ and 
K2 are the solutions of the characteristic equation 

an 2 + (2biiu + 2/3 x s)k + 2i(3 v ujs - s 2 - cw 2 = 0. (3.43) 

Without restriction we can assume a = —1. Moreover, since the sign of !Rk does not depend on uj, we set w = 0. 
We then obtain 



K h2 =p x s±^(f3 x2 -l)s 2 



For !Rs > 0, we have !Rki,2 > and there is no bounded solution u. Therefore no boundary condition is needed and 
the problem is well-posed. □ 

As difference approximation for the half-plane problem, we can use either l|3.9|l or 13.30JI . We study the stability of 
the approximations by mode analysis. Below we show that 113. 9J1 is stable with pi(D) in H3.12I) . The stability of the 
other approximations with p\ (D) in IJ3.11H and IJ3.13J1 can be shown in the same way. 

On a uniform spatial grid flh = (yh, fih), v = 0, 1, 2, . . . , [i, = 1, 2, . . . , N, with spacing h, let v(t) := u vll {t) be the 
gridfunction approximating u(x u , y M ,£). We consider the shifted wave equation H3.36I) and approximate it by 

v u = 2{(3 X D +X + (3 y D +y )v t - {I3 X D +X + [3 y D +v ) 2 v + ( ai D +x D_ x + 2b 1 D 0x D 0y + Cl D +y D_ y )v, (3.44) 

for v = 1, 2, ... . For every fixed /i, we need one extra boundary condition to determine uq^. We use a third order 
extrapolation 

h 3 D 3 + x u 0fi = 0. (3.45) 

We consider bounded solutions of type 

«„„(*) = e st +^V„ IMU<oo. (3.46) 

Putting IJ3.46II into IJ3.44II . we get the eigenvalue problem 

ip v s 2 - 2—(Lp v+1 - ip v )s - 2— (isin£ - 2 sin 2 -)tp v s 

8 x2 8 x B y £ B y2 £ 

+ -r^isPv+i - 2^+1 +<p I/ ) +2 (<fv+i - <p„)(ism£- 2 sin 2 -) + — ^-(isin 2 ^- 2 sin 2 -) 2 <p v 

~ T2&V+ 1 ~ 2< ^ + <P»-i) - p-isin^^+i - (fiu-i) +4-p- sin 2 -tp v = 0, £ = u>h. (3.47) 

The approximation J3.440 - I3.45t is stable if and only if the Kreiss condition is satisfied, or equivalently if l|3.47f) has 
no eigenvalue s with !Rs > J2lJ • The constant-coefficient ordinary difference equation 113. 4711 has solution of the form 

3 

j=i 
where Kj are the three solutions of the characteristic equation 



- 2( t —{K-l) + t —{ism£,-2sin 2 ^)\s+ I ^-(k - 1) + ^-(isin£ - 2siii- ,, 

ai (k — l) 2 6i , 1 . ci o £ , 

^L_^__ (K __ )iSlne + 4 _sm 2 |^0. (3.48) 
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By Lemma 12.1.6 of J2J|, for 5fts > the characteristic equation i!3.48|l has no solutions with \k\ = 1 and there is 
exactly one solution with |k| < 1. Roughly speaking, the number of left points in the difference stencil determines 
the number of solutions to the characteristic equation with |«| < 1. We call this solution K\ and write the bounded 
solution as 

Uv (t) = e st+ ^Vi/<. (3.49) 

By substituting l|3.49|) into the boundary condition Ij3.45f) . we get 

a 1 (Ki-lfe st+iu "* h =0. (3.50) 

Since m ^ 1 for !Rs > 0, 113. 50|) has only the trivial solution o\ — 0. Now, we let k — ► 1 and investigate if there is any 
sequence {s} such that 5Rs — ^ with 5Rs > 0. We then get from H3.48JI 



x2 r>ay/--c o„;„2^\- , o»2/- • <■ i„;„2£n2 , , „ ■ 2 £ 



2/3 a (isin£ - 2 sin 2 -)s + /J y2 (isin£ - 2 sin 2 -) 2 + 4c x sin 2 - = 0, s = sh, (3.51) 



and therefore 



.2 £\, . / „„ -2 £ 



§ = /3 y (isin£-2sm 2 p± W-4c x sin 2 |. (3.52) 

Since /3 y > and Ci > 0, we have 5Rs < if £ -y> 0. In the case where £ — ► 0, we get from 113.4811 



s 



.^.. 1) + ^ (/ ,. ,,._*. («-)• 



2^- s ( K -l) + ^— ( K -l) 2 --^ i- = 0. (3.53) 



Letting s — -> 0, we then get from l|3.53|l that K\ t 2 = 1 and K3 = ai//3 x < 1. Since for 5Rs > there is no solution with 
|k| = 1, the only solution is K3 which is strictly less than 1 and does not converge to 1. Therefore there is no positive 
sequence {s} such that 3?s — > for |£| < n. Now, we can prove the following theorem: 

Theorem 6. The approximation H3.44jl - lj3.45jl is stable. 

Proof. Since there is no eigenvalue s with 5Rs > to the eigenvalue problem H3.47JI giving bounded solutions l|3.4fill . 

the Kreiss condition is satisfied and stability follows. □ 

IV. TESTS OF THE SUPERLUMINAL ALGORITHMS 

In the subluminal case where the evolution proceeds in a timelike direction, the W^-algorithm jl.8|l provides an 
accurate, flux-conservative, second order treatment of the IBVP. This was proved for a ID quasilinear wave equation 
in [lj| using the discrete energy method. In [lfiL Il2l|. the results were extended to the 3D case and applied to the 
harmonic Einstein system ill,2fl . The semi-discrete conservation laws extend to the principle part of the harmonic 
Einstein system and contribute to excellent long term performance in test problems. We use this Vy-algorithm to 
treat the outer region of the model excision problem considered in Sec. 

In this model problem, the inner boundary is chosen to be spacelike, corresponding to the strategy for excising an 
interior singularity. The evolution near the inner boundary proceeds in a spacelike direction (superluminal shift) so 
that the spatial grid tracks the boundary. For this superluminal case, the W-algorithm is unstable and one of the 
algorithms considered in Sec. IIIII must be used. These algorithms are either given by l|3.9|l . with Pi{D) given by one 
of the approximations H3.11j) - (I3.13II . or by H3.30II . 

In the case of the 2D shifted wave equation Ijl.fijl , the choice i|3.11JI reduces to the centered algorithm 



V := f (ft - /3 x D 0x - l3 v D 0y ) 2 - cnD +x D_ x - aD +y D_ v - 2b x D Ox D 0y )u = 0; 
the choice i|3.12H reduces to 

V+ := ( (fi t - X D +X - P y D +y f - ai D +x D_ x - dD +y D_ y - 2b 1 D 0x D 0y )« = 0, (4.2) 



in which the shift terms are treated by first order accurate one-sided difference operators; the choice IJ3.131 reduces 
to 



V p := (d t - (3 x D px - P v D py ) z - ai D +x D_ x - c x D +y D_ y - 2b 1 D 0x D y )u = 0, (4.3) 
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in which the shift terms are treated by second order accurate one-sided difference operators IJ3.14H : and 113,3011 is 
related to the subluminal VF-algorithm 11.8JI by 



V a := W + !L ( ai (D +x D- x ) 2 + a 2 (D +y D- y )Au = 0, (4.4) 

where Theorem 4 guarantees stability provided the inequalities 

2 

a\>j3 x —ai = —a, a± > 0, (4.5) 

u 2 >P y2 -ci =-c, a 2 >0, (4.6) 

are satisfied. 

The T^-algorithm is related to the VF-algorithm by the second order accurate modification 

V = W + !j (V (D +X D_ X ) 2 + py 2 {D+yD-yA u = 0. (4.7) 

In the subluminal case where the W and V algorithms can be compared, tests show that the 14^-algorithm has 
considerably better accuracy due to its more compact stencil [l{|]. Here we carry out a set of 2D superluminal tests 
to compare the performance of the superluminal algorithms in a periodic test problem (smooth toroidal boundary 
conditions) where the effect of the boundary is eliminated. The first order accurate V+-algorithm 14.21 is highly 
dissipative and much less accurate than the second order accurate V p version ll4..'"it . For these reasons, we restrict our 
test comparisons to the V, V p and V a algorithms. 

2 2 

The F-algorithm is a special case of the V Q -algorithm (I4.4f) where a\ = (3 X and a 2 = (i v . The accuracy of 
the V^-algorithm might be expected to depend on the relative weight of the higher order terms responsible for the 
stretched stencil in 14.411 . For example, the V^-algorithm might be expected to be most accurate for the minimum 
values, ai = (M — a )/2 and a 2 = (|c| — c)/2, which are allowed by 14. -'it and 14.611 . This is true for the case when a 
and c are positive, for which a\ = a 2 = and the V^ -algorithm reduces to the W-algorithm. 

When a < is negative and c > 0, the optimal value for a 2 remains but the optimal value for a\ is not necessarily 
the minimum allowed value a\ = —a. This value would result in approximating the ad x term in the wave operator by 
o.Dq x , which decouples the even and odd grid points. Although the optimal choice of a\ in this case is not obvious, 
the combination a\ = (3 X and a 2 — would give better accuracy than the ^-algorithm. No general guidelines are 
suggested by examining the truncation error in the V^-algorithm, which to order h 2 is given by 

r = ^ f(a - 3«i)^ + (c - 3a 2 )a* + ±b(d 3 x d y + d x d a y ) + 4(3 x d 3 x d t + 4(3 x d%d t ) u. (4.8) 

Note that the values a\ = a/3 and a 2 = c/3 correspond to the fourth order accurate approximations to the terms 
ad 2 and cd 2 in the wave operator. However, these choices are not allowed in the superluminal regime, where stability 
requires a, > 0. 

As a test problem for comparing the accuracy of these evolution algorithms in the superluminal regime we pick a 
case where both a and c are negative. We consider the wave equation 

d 2 + A(d x + d y )d t ~ 3d 2 ~ 3d 2 - 8d x dy] u = 0. (4.9) 



which arises from a 2D version of 112.411 with shift (3 X = [3 y = 2. With this superluminal choice of shift, there are no 
characteristics in the (x > 0,y > 0) directions. Waves propagating along the diagonal have the form 

u = F[x + y + (4 + y/2)t] + G[x + y + (4 - V2>]. (4.10) 

In our test, we simulate the solution 



U = sin ( 2Tr[a; + y + (4 + V2)t] j 



(4.11) 



in the domain —.5 < (x,y) < .5, on a grid with N = 
particular solution, the symmetries d x u = d y u = dtu/(A 
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200 points, with periodic boundary conditions. For this 
- \/2) imply that the truncation error l|4.8|l has a minimum 



at a\ 



tti 



a m , where 



13 + 8^ 



.1045695. 



(4.12) 



FigureQlplots the £oo norm of the numerical error in the scalar field obtained in the simulation of (14. 1 If) by evolving 
the wave equation 114. 9J1 with the Va-algorithm, for various values of a\ = «2 = ot. The error for a = 8.1045695 is 
extremely small and the plots confirm that a m is indeed the optimal value. The value a = 4 corresponds to the 
^-algorithm, which gives significantly larger error. The value a = 3, which is the smallest value allowed by stability, 
gives even larger error. 
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PIG. 1: The loo norm of the error of the scalar field obtained with V a algorithm on a grid of 200 points is plotted vs time, in 
the interval < t < 1. For the value a m ~ 8.1045695, the error is barely discernible and the plots clearly indicate that a m is 
the optimal value. The value a = 4 corresponds to the ^-algorithm, which has significantly larger error. 

The error in Fig. Q] is predominantly phase error. Figure [2] shows snapshots of u(t = 100, x) (100 crossing times) 
for the simulation of (|4.11l) using the V, V p and V a algorithms, with a = 8. The simulations are compared with the 
analytical solution at t = 100. The solution with the V p - algorithm leads in phase while that with the ^-algorithm 
lags in phase and it has slightly better accuracy. As expected from the above error analysis, the V^-algorithm, with 
a = 8, is extremely accurate. 
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FIG. 2: Snapshots of the scalar field u(t = 100, x) obtained with the V, V p and V a algorithms, compared with the analytic 



solution. The phase error with the Vp- a lg°rithni i s larger than with the V-algorithm. The Va- algorithm, for a ■ 
accurate and barely distinguishable from the analytic solution. 



8, is extremely 



SIMULATION OF A MODEL 2D EXCISION PROBLEM 



In this section, we simulate a simple 2D model of the excision problem in which the inner boundary S is spacelike 
and the outer boundary T is timelike, with a horizon ri in between. In the inner region between S and ri, since the 
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shift is superluminal, the operator P in i|3,23|l is non-elliptic and both characteristics leave the inner boundary. In the 
outer region between H and T, since the shift is subluminal, the operator P is elliptic and one characteristic leaves 
T and the other enters T. 

To model a wave pulse propagating into a horizon, we consider the shifted wave equation with a source term F, 



u t t = 2{f3 x u xt + P v u y t) + au xx + 2bu xy + cu yy + F(x, y, t), 



(5.1) 



on the spatial domain (x, y) £0= [-2,2] x [-2,2], and t > 0. We set the coefficients (3 X — (3 V — 2, a = 0.5(x — sin^), 
b = 0.5 and c = 5, for which the problem is well-posed. The spacelike boundary S at x = —2, the timelike boundary 
T at x = 2 and the horizon H are shown in Fig. 03 The horizon satisfies ac — b 2 = 0, which determines the curve 



x= 0.1 + sin—. 



(5.2) 




T 



FIG. 3: Computational domain with the spacelike boundary S on the left, the timelike boundary T on the right and the 
sinusoidal shaped horizon TL in between. The solution is periodic in the vertical y-direction. 



For the smooth function 
2 



F(x, y,t) = -^(-(l-2/3 x - a)(u -2{t + x- x f) - 4(/3 a + b)(t + x- x )y + c(a - 2j, 2 )) e 
the equation Ij5.1|l has the solution 



u(x,y,t) = e 



(H--»o) / +» / 



(5.3) 



(5.4) 



which is a left-traveling wave packet, initially centered about x = 0.5 outside the horizon and propagating towards 
the spacelike boundary. Here we set a = 0.05. We uniformly discretize the spatial domain as x v = vh and y^ = fih 
with v, fi — 0, ±1, . . . , ±N with the grid size h = j?. 

The global simulation of the model problem in the region between S and T is carried out by combining the 
superluminal ^-algorithms established in Sec. IIIII with the subluminal IF-algorithm. The spacelike boundary and the 
superluminal region are treated with one of the ^-algorithms. A region containing the timelike boundary is treated 
by the VF-algorithm. 

We consider the following three global algorithms: 

• Algorithm 1. The superluminal region is treated by the ^-algorithm 14.3JI . In the subluminal region we use 
the VF-algorithm. We introduce a cut-off function </> which is when a > and c > and is 1 when a < or 
c < 0. Then we use the following approximation 

cj)V + {I - <j))W = Q. 

• Algorithm 2. This is similar to 1, except the superluminal region is treated by the VJ,-algorithm IJ4.3II . which 
is then blended to the W-algorithm in the same way as in 1. 

• Algorithm 3. We use the V^-algorithm J4.4II with a± = (\a\ — a)/2 and «2 = (|c| — c)/2. 
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The initial data and boundary condition at x = 2 are chosen according to the exact solution i|5,4|l . In the first 
and third algorithms, we need two extra boundary conditions at v = —N,—N + 1. In the second algorithm, we 
need only one boundary condition at v = —N. We use third order extrapolations as the extra boundary conditions. 
In the y-direction we use periodic boundary conditions. For the integration in time, we use the standard 4th order 
Runge-Kutta method. 

Figure Bl shows the initial wave pulse and the pulse at a later time t = 2 computed by the third algorithm. 
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(a) The initial wave pulse at t = 0. 



(b) The wave pulse at t = 2. 



FIG. 4: A pulse propagating across the horizon. The left figure shows the initial pulse in the region outside the horizon. The 
right figure shows the pulse at a later time, after it has crossed the horizon and is incident on the inner spacelike boundary. 



For the gridfunction u„ M (t) approximating u(x^,y ll ,t), we define the discrete norm as 

N 
v,)j,=—N 

where h = Ax = Ay is the gridlength. We then define the convergence factor by 

C(t) = log 2 ( iS^J . *(*) = uC^Wm.*) " «vm(*)> 

where £{t) is the error at time t, and u(xi/,j/^,t) is the exact solution computed by H5.4H . 

Figure shows the norm of the error versus time for the three algorithms with h — 0.02 and At = 0.001. 



(5.5) 



(5.6) 
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FIG. 5: Norm of the error ||£|k versus time, with h = 0.02. 



Figure shows the convergence factor as a function of time for the three algorithms with h = 0.04 and At = 0.001. 
It confirms the second order accuracy of the algorithms in space. The jumps in the convergence factor at about t = 2 
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is a result of using third order extrapolations at the spacelike inner boundary, while we use second order evolution 
algorithms. At this time the pulse reaches the spacelike boundary and an increase in the order of accuracy, from 2 to 
3, is expected. 




(a) Algorithm 1. 



(b) Algorithm 2. 



(c) Algorithm 3. 



FIG. 6: Convergence factor C versus time indicates that the algorithms are second order accurate in space. 

The third algorithm gives a better accuracy than the other two. The second algorithm, in which we use the second 
order one-sided stencil with extrapolation in one ghost point, gives a slightly smaller error than the first algorithm, 
in which the second order centered stencil with extrapolation in two ghost points is used. 
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